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Abstract 

Ratio-dependent predator-prey models have been increasingly favored by field ecologists where 
predator-prey interactions have to be taken into account the process of predation search. In this 
paper we study the conditions of the existence and stability properties of the equilibrium solutions 
in a reaction-diffusion model in which predator mortality is neither a constant nor an unbounded 
function, but it is increasing with the predator abundance. We show that analytically at a certain 
critical value a diffusion driven (Turing type) instability occurs, i.e. the stationary solution stays 
stable with respect to the kinetic system (the system without diffusion). We also show that the 
stationary solution becomes unstable with respect to the system with diffusion and that Turing 
bifurcation takes place: a spatially non-homogenous (non-constant) solution (structure or pattern) 
arises. A numerical scheme that preserve the positivity of the numerical solutions and the bound- 
edness of prey solution will be presented. Numerical examples are also included. 

Key words: reaction-diffusion system, population dynamics, bifurcation, pattern formation. 
PACS: 35K57, 92B25, 93D20. 



1. Introduction 

Since it is rare to find a pair of biological species in nature which meet precise prey-dependence 
or ratio-dependence functional responses in predator-prey models, especially when predators have 
to search for food (and therefore, have to share or compete for food), a more suitable general 
predator-prey theory should be based on the so-called ratio-dependent theory (see [H, 0, 0, 0] ) • The 
theory may be stated as follows: the per capita predator growth rate should be a function of the 
ratio of prey to predator abundance, and so should be the so-called predator functional response. 
Such cases are strongly supported by numerous field and laboratory experiments and observations 
(see, for instance, [Ijg, Q, 

Denote by N(t) and P(t) the population densities of prey and predator at time t, respec- 
tively. Then the ratio-dependent type predator-prey model with Michaelis-Menten type functional 
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response is given as follows: 
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~dt 
dP 
~dt 
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Q(P) + 



aNP 
mP + N' 
bN 



mP + N 



(1.1a) 
(1.1b) 



where a,b,m,K, and r are positive constants. In (jl.l|) . Q(P) denotes a mortality function of 
predator, and r and K the prey growth rate with intrinsic growth rate and carrying capacity in 
the absence of predation, respectively, while a, b, and m are model-dependent constants. 

From a formal point of view, this model looks very similar to the well-known Michaelis-Menten- 
Holling predator-prey model: 



dN 
~~dt 
dP 
~dt 



rN 
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Q(P) + 



aNP 
~ c + N' 
bN 



c + N 



(1.2a) 
(1.2b) 



Indeed, the only difference between Models (jl.ip and (|1.2p is that the parameter c in (jl.2p is 
replaced by mP in (jl.ip . Both terms mP and c are proportional to the so-called searching time 
of the predator, namely, the time spent by each predator to find one prey. Thus, in the Michaelis- 
Menten-Holling model (|1.2p the searching time is assumed to be independent of predator density, 
while in the ratio-dependent Michaelis-Menten type model (jl.ip it is proportional to predator 
density (i.e., other predators strongly interfere). 

Predators and preys are usually abundant in space with different densities at difference posi- 
tions and they are diffusive. Several papers have focused on the effect of diffusion which plays 
a crucial role in permanence and stability of population (see 0, Eol . 11. 12. QT 



references therein). Especially in 
was extensively studied, and in [11 
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14, UM, and the 



the effect of variable dispersion rates on Turing instability 
the dynamics of ratio-dependent system has been analyzed in 
details with diffusion and delay terms included. Cavani and Farkas (see (l6| ) have considered a 
modification of (jl.2p when a diffusion was introduced, yielding: 



dN A7 / N 

— = r N 1 

dt V K 



aNP „ d 2 N 



dP 
~dt 



P 



-Q(P) + 



c + N 
bN 



c + N 



+ D. 



dx 2 

d 2 p 



dx 2 



x G (0,l),t > 0, 



where the specific mortality of the predator is given by 

J + SP 



Q{P) 



1 + P ' 



(1.3a) 
(1.3b) 

(1.4) 



which depends on the quantity of predator. Here, the positive constants 7 and 5 denote the 
minimal mortality and the limiting mortality of the predator, respectively. Throughout the paper, 
the following natural condition 

0<7<5 (1.5) 
will be assumed, and we will consider the case of the constant diffusivity, Dj > 0, i = 1,2. The 
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advantage of this model is that the predator mortality is neither a constant nor an unbounded 
function, but still it is increasing with the predator abundance. On the other hand, combining 
(jl.ip and (jl.3p , many authors (see 17J, [l5|, [18|] , for instance) have studied a more general model as 
follows: 



8N 
~dt 

dP_ 

dt 



rN 1 



P 



N 
K 

-Q(P) + 



aNP 
mP + N 
bN 



d 2 N 

fDi— , x E (0,1), t > 0, 



+ D. 



dx 2 

d 2 p 



dx 2 ' 



x G (0,l),t > 0, 



mP + N 

with the specific mortality of the predator somewhat restricted in the form 

Q(P) = d. 



(1.6) 



In this paper we consider a ratio-dependent reaction-diffusion predator-prey model with Michaelis- 
Menten type functional response and the specific mortality of the predator given by (jl.4p instead 
of (jl.6p . We study the effect of the diffusion on the stability of the stationary solutions. Also 
we explore under which parameter values Turing instability can occur giving rise to non-uniform 
stationary solutions satisfying the following equations: 



dN 
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dP_ 

dt 
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aNP 
mP + N 
bN 



l + P mP + N 



+ Dy 

+ D 2 



d 2 N 
dx 2 ' 
d 2 P 

dx 2 ' 



x G (0,1), t > 0, 
x G (0,l),t > 0, 



(1.7a) 
(1.7b) 



assuming that prey and predator are diffusing according to Fick's law in the interval x G [0,/]. We 
are interested in the solutions N, P : (1,0) X M + — > M + fulfilling the Neumann boundary conditions 



N x (0,t) = N x (l,t) = P x (0,t) = P x (l,t) = 0, 



li 



and initial conditions 



N(x, 0) > 0, P(x, 0) > 0, x G (0, 1). 
For simplicity, we nondimensionalize the system (jl.7|) with the following scaling 



~ N 
t = rt, N = -, 



P 



mP 



and letting 



a 



a „ 7 ~ S b K D x , D 2 

' 7 = t> " = Ti e= -, /3 = —, dx = — , d 2 = — • 

mr b o r m r r 



For the sake of simplification of notations, dropping tildes, the system (jl.7p takes the form 



dt K ' P + N 



aNP , d 2 N . JN 

+ d 1 ^ r , x G (0,1), t > 0, 



dx z 
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7 + S/3P 

' 1 + (3P 'P + N 
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8 2 P 
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(1.9a) 
(1.9b) 
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Set 



where 



F 2 



u := 



N 
P 



D :-- 



di 
d 2 



Fi(JV,P) = N(l - N) 



aNP 
P + N' 



F 2 (N,P)=eP 



7 + 60P N 
' 1 + PP + P + N 



Then the system (II. 9|) with the boundary conditions (jl.8p takes the form 



<9 2 u 

u t = F(u)+D^; u^O,*) =u x (/,t) = 0. 



;i.io) 



Clearly, in case the predator and prey are spatially homogeneous, the spatially constant solution 
u(t) = (N(t), P{t)) T of (jl.lOp . fulfilling the boundary conditions obviously, satisfies the kinetic 
system 

u t = F(u). (1.11) 



2. The model without diffusion 

In this section we will study the system fll.9)l without diffusion, i.e. 

aNP 



dN 
~dt 



N(l - N) 



P + N' 



d 4 = ^ 

dt 



+ 
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1 + 5I3P 

' 1 + PP 'P + N 



(2.1a) 
(2.1b) 



In particular, we will focus on the existence of equilibria and their local stability. This information 
will be crucial in the next section where we study the effect of the diffusion parameters on the 
stability of the steady states. 

The equilibria of the system (12. ip are given by the solution of the following equations 



N(l - N) 



aNP 
P + N 



0, eP 



+ 



N 



' 1 + PP 'P + N 



0. 



The system has at least one equilibrium with positive values. This is the point of intersection of 
the prey null-cline 

(1 - N)N 



P = H l (N) 



a-(l-JV) 



and the predator null-cline 



P = H 2 {N) 



7 - P(l - 6)N + 2^/{ 7 - p{l - 5)N} 2 - 4/35(1 - j)N 
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Thus, denoting the coordinates of a positive equilibrium by (N,P), these coordinates satisfy P 

H 1 (N)=H 2 (N). 

The Jacobian matrix of the system (|2,ip linearized at (N, P) is 



,4 



6i -0 2 

03 "04 



(2.2) 
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where 
and 



trace A = Q 1 - 4 , det A = 2 3 - ©i© 4 



aNP „ aN 2 
0i = -N + -=-=-, 2 



(P + iV) 2 (P + jV) 2 

eP 2 e0P(S—r) , eiVP 

^3 = 7=^ ^T7> fc) 4 = — ^=77: + 



(P + N) 2 (1 + /3P) 2 (P + N) 2 

The characteristic equation is given by 

A 2 - (trace A) A + det A = 0. 

Recall that (N, P) is locally asymptotically stable if Re A < 0, which is equivalent to have trace A < 
and detj4 > 0. For this, we will assume that 

01 < 4 , 02©3 > 0104- (2.3) 
Remark 2.1. Due to (|1.5|) . we see that 4 > 0. If Qi < 0, then the two conditions in (|2.3|) hold. 



3. The model with diffusion 

In this section we will investigate in Turing instability and bifurcation for our model problem. 
We will also study pattern formation of the predator-prey solutions. 

3.1. Local existence of solutions 

Before studying the stability of equilibrium solutions, we will discuss about the local existence 
and uniqueness of solution for a given ratio-dependent reaction-diffusion predator-prey model. 
Applying the criteria for the local existence of solution (see [jjj, l2fj| ) to the nonlinear parabolic 
systems (jl.lOj) . we see that there exists a unique local solution of the given system. 

Let Q. be a bounded region in M. n ,n > 2, with smooth boundary d£l and v denotes the unit 
outward normal to Q. Then Morgan considered in reference ([lii]) essentially of the form 

u t (x,t) = DAu(x,t) + /(u(x,i)), x£n,t>0, (3.1a) 

du( f>^ = o, xedn,t>o, (3.ib) 

av 

u(x, 0) = u (x), (3.1c) 

where u:!lx(0, oo) — > M m , / : M m — > R m is a locally Lipschitz continuous function, D is an m x m 
diagonal matrix with diagonal entries dj > 0, and uo : O — > M m is bounded and measurable. Then 
the following theorem holds (li| : 



Theorem 3.1. Under the assumptions on (13. If) stated above, there exists T max > and M = 
(Mf) e C([0,T max ),IR m ) suc/i toot 

(i) (E2]) /ios a unique classical solution u on!lx [0, T m ax) which cannot be continued to [0, T) /or 
any T > T max , and 

(ii) lujC-.tJIocn < Mj(t) for all 1 < j < m,0 < * < T max . 

Moreover, i/T max < oo, i/ien |uj(-, i)|oo,c — > oo as i — > P max - /or some 1 < j < m. 
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Defining iq(0, 0) = and ^(0,0) = in our model (jl.lOp . Theorem 13.11 implies local existence 
and uniqueness. More precisely, there exists T max > and Nm and Pm £ C([0,T max )) such that 

(i) (|3.ip has a unique classical solution u = (N, P) T on [0, I] x [0, T max ) which cannot be continued 
to [0, T) for any T > T max , and 

(ii) lA^^L^o,/) < N M (t) and |P(-, i)^,^) < P M (t) for < t < T max . 

Moreover, if T max < oo, then either \N(-, i)|oo,(o,i) ^ oo or |P(-, t)|oo,(o,0 -> oo as i -> T max -. 



5.5. Turing instability 

Definition 3.2. VFe say that the equilibrium (N,P) is Turing unstable if it is an asymptotically 
stable equilibrium of the kinetic system 112.1)) but is unstable with respect to solutions of M.9)) (see 

0/;. 

An equilibrium is Turing unstable means that there are solutions of (jl.lOp that have initial 
values u(x, 0) arbitrarily closed to u (in the supremum norm) but do not tend to u as t tends to 
oo. 

We linearize system (jl.9p at the point (N,P): setting v = (v\,V2) T = (N — N,P — P) T , the 
linearized system assumes the form 

dx 1 

while the boundary conditions remain unchanged: 



vt = Av + D—^, (3.2) 



v*(0,t)=v x (Z,t) = 0. (3.3) 

The linear boundary value problem (|3.2p - (|3.3p can be solved in several ways. In particular, the 
Fourier's method of separation of variables assumes that solutions can be represented in the form 
v(x, t) = i/t(x)y{t), with y : [0, oo) M 2 , if) : [0, 1] -> R. Then 

d i = (A-(D)y, (3.4) 

and 

-Vxx = CV', ^*(0) = V*(0 = 0. (3.5) 
The eigenvalues of the boundary value problem (|3.5p are 

V 



= ( J ~T ) » J = 0,l,2,--- (3.6) 



ipj(x) = cos — — . (3-7) 



with corresponding eigenfunctions 

J7TX 
1 

Clearly, = Co < Ci < C2 < ' • • • These eigenvalues are to be substituted into (|3.4|) . Denoting by 
yij and y2j the two linearly independent solutions of (|3.7|) associated with C = Cj> the solution of 
the boundary value problem (|3.2p - (|3.3p is obtained in the form 

00 

V (M) = ^(aij-yij(*) + 02jy2i(t))cos^ (3.8) 
3=1 
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where aij, i = 1, 2, j = 0, 1, 2, • ■ ■ , is to be determined according to the initial condition v(x, 0). For 
instance, if yii (0) = (1, 0) T , y 2j (0) = (0,1) T for j = 0, 1, 2, • • • , 



aio 

O20 



v(x, 0)dx, 



&2k 



A; = 0,1,2, 



Set 



2 f kirx 
- / v(x, 0) cos — — ax 

o 

B(C)=A-CA Bj = B(Q) = A — QD. 



(3.9) 



According to Casten and Holland 10(j, if both eigenvalues of Bj have negative real parts for all j, 
then the equilibrium (iV, P) of (jl.lOj) is asymptotically stable; if at least one eigenvalue of a matrix 
Bj has positive real part, then (N,P) is unstable. Recalling (|2.2p . the trace and determinant are 
given by 



trace Bj 
det Bi 



e 1 -e 4 -Ci(rfi + d 2 ). 

G 2 6 3 - 9i9 4 + W©4 - d 2 6i} + C|did 2 . 



(3.10a) 
(3.10b) 



Notice that f|2.3[) implies that trace < 0. Therefore the eigenvalues of Bj have negative real parts 
if det Bj > which is guaranteed in case 

di9 4 > d 2 &i, (di0 4 - ^ 2 ei) 2 - 4^1^(0203 - 0i0 4 ) < 0. 

Notice that det-B^ < for all sufficiently large j if di© 4 — d 2 0i < 0, since the eigenvalues Q is 
monotonic increasing with its limit oo. Therefore, one has the following theorem: 



Theorem 3.3. Assume that (|1.5p and (|2.3p . Then the equilibrium point (N,P) of (jl.lOp is 

asymptotically stable if 



di© 4 > d 2 0i, (di6 4 - rf 2 ©i) 2 - 4^^(0203 - 0!0 4 ) < 0; 
while it is Turing unstable if 

di6 4 - d 2 ©i < and (cZi0 4 - d 2 ©i) 2 - 4did 2 (0 2 3 - ©i© 4 ) > 0, 
or if there exist a positive integer k such that 

detB k = 2 ©3-0i0 4 + CfcR04-d20i} + C^i^2<O. 



(3.11) 



(3.12) 



(3.13) 



3. 3. Pattern formation 

For a nonnegative real parameter A consider the reaction-diffusion system to find u : (0, 1) x 
(0, oo) -» M n such that 

u t = F(u;A) + £>(A)|^, (3.14) 
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where D is a non-negative diagonal matrix depending smoothly on A and F : M n x [0, oo) — » W 1 is 
a smooth function. Suppose (|3.14p is equipped with the Neumann boundary condition 

u x (0,t) = u x (l,t) =0. (3.15) 

Assume further that for some u € W 1 we have F(u; A) = for all A € [0, oo), i.e. u is a parameter- 
independent constant stationary solution of (|3.14p - (|3.15l) . 

Definition 3.4. We say that u undergoes a Turing bifurcation at Ao £ [0, oo) if the solution u is 
asymptotically stable for < A < Ao, while it is unstable for Ao < A, (or vice versa, i.e. the regions 
for asymptotical stability and instability may be exchanged), and in some neighborhood of Ao the 
problem \3. ll$ - $3. 13\) has non-constant stationary solution (i.e. solution which does not depend on 
time but depends on space.) 

With d\ fixed, regarding d 2 as the parameter A, we will consider the linearized system (j3. 2|) - 
(|3.3[) as a parameter-dependent problem in the setting (|3.14|) - (|3.15|) . Notice that u(x,t) = (0, 0) T 
is clearly a solution for (|3.2p - (|3.3p . Then the condition for a Turing bifurcation for the linearized 
system (|3.2p - (|3.3|) is given as follows: 

Theorem 3.5. Suppose that traced < and det A > 0. 

di > % (3.16) 



then the zero solution of the linear problem 113. ty) - [373\) is asymptotically stable for all d2 > 0. 
(ii) U 

^<di<% (3.17) 
then the zero solution of the linear problem ^3.2) - ^3~3) undergoes a Turing bifurcation at 

, , e 2 e 3 -e 1 e 4 + Cirfie 4 ,„ 1fi . 

«2 :- a 2c rit -—Tf\ TT~\ 1 (3. 18 J 

Proof, (i) Rewriting (|3.10b[) as 

det Bj = 9 2 e 3 - 6i8 4 + 0^i@4 - (jd 2 (@i - Cjdi), 

we see from (ll.5p and (I2.3P that det Bj > for all j = 0, 1,2, ••• if d\ > ©i/Ci holds, since 
CjiJ = 0, 1,2, ••• forms a monotone increasing sequence (|3.16p . Therefore, the zero solution of 
(I3.2|) - (l3.3p is asymptotically stable under such conditions. 

(ii) Suppose d\ satisfies (13.17j) and choose A = d 2 as given in (I3.18|) . Then deti?i = 0. Clearly, 
we have deti?i > for < d 2 < d 2cr i t , and deti?i < for d 2cr i t < d 2 . In both cases det-Bj > 
0,j 1. Again by Casten and Holland [k| as quoted just after formula (|3,9p . the zero solution 
is asymptotically stable for < d 2 < d 2cr it, and it is unstable for d 2cr u < d 2 . If d 2 = d 2cr n, one 
eigenvalues of B\ becomes zero and the other is trace Si, which is negative. Denote the eigenvector 
corresponding to the zero eigenvalue by yn = (r]i,r] 2 ) T , i.e. 

Biyn = {A- Ci£>)yn = 0, yn + 0. 
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As we can see from (|3.4p - (|3.7p the function 



vi(x,t) := ynV'iO) 



m 
m 



TTX 

cos- 



is a spatially non-constant stationary solution of the linearized problem (|3.2p - (|3.3p . This implies 
that the zero solution undergoes Turing bifurcation at efecrit- This completes the proof. □ 

In the remaining part of this section we will extend the latter result about the Turing bifurcation 
of the zero solution of the linearized system to the non- linear problem (|1.10p . For this we need the 
following: 

Theorem 3.6. Let X and Y be Banach spaces, U = V x S an open subset of X x R, and 
f G C 2 (U; Y)such that f (0, A) = 0, A G S C R. Denote by L w = f v (0, A ) and L 12 = f v ,A(0, A ) the 
linear operators obtained by differentiating f with respect to its first variable only and the first and 
second variables at v = G V, Ao G S, respectively. Assume that the following conditions hold: 

(i) the kernel of L\q, the subspace N(L\q) of X is a one- dimensional vector space spanned by 
vi G X; 

(ii) the range of L\q, the subspace R(L\q) ofY has codimension 1, i.e. dim[Y '/ 'R(L\q)\ = 1; 

(iii) L 12 vi i R(L W ). 

Let Z be an arbitrary closed subspace of X such that X = [Span v±] © Z; then there is a 5 > 
and C x -curve 0, A) : (-5, 5) -> Z x S such that; (f)(0) = 0; A(0) = A ; f (svi + 8(/>(s), A(s)) = for 
\s\ < 5. Furthermore, there is a neighborhood of (0, Ao) such that any zero of f either lies on this 
curve or is of the form (0, Ao). 

Proof. The idea of the proof is to introduce a new parameter s which enables to apply immediately 
the implicit function theorem for the function F G C 1 (U x Z,Y) defined by 



F(A,a,*) 



'±f(svi + sz,A) if s ^0, 
Lio(vi+z), if s = 0. 



See, for the details of the proof of the theorem, pp. 172-173 of 2l[. □ 

Remark 3.7. In what follows the role of the space X will be played by 

X = {VGC 2 ([0,/];R 2 ) : V E (0) = V x (l) = 0} (3.19) 

with the norm \\f\\x = ]Co<a<2 su Pxg[o,/] | c^^f (a?) | , where \ ■ \ denotes the usual vector or matrix- 
norm, while Y = C°([0,/],R 2 ) with the norm ||f||y = sup^g^^] |f(a:)|- However, in choosing the 
subspace Z of X we shall use the orthogonality induced by the inner product 

i 

(v, w) = / [vi(x)wi(x) + V2[x)w2{x)\ dx, for v = (v±, V2) T , w = (wi, W2) T ■ 



Theorem 3.8. Suppose that traced < and detvl > 0. 

(i) If A3. 1 6\) holds, then the constant solution u = (N, P) T of the nonlinear problem il.lO\) is 
asymptotically stable. 



9 



(ii) If (0, r]2) T is not parallel to the second eigenvector y 2 i of B\ and d\ satisfies \3.11 ), then at 
di = dicrit the constant solution u undergoes a Turing bifurcation. 



Proof, (i) follows immediately from the asymptotic stability of the zero solution of the linear prob- 
lem (pT2D - (nOj) . 

(ii) As in the proof of (i) of Theorem 13.51 we have that u is asymptotically stable for d% € 
(0,d2crit), while it is unstable for <i 2 E (d 2cr jt, oo). We have to show the existence of a stationary 
non-constant solution in some neighborhood of the critical value d^crit- Such a stationary solution 
u satisfies the following system of second-order partial differential equations 

Du xx + F(u) = 0, u x (0) = u x {l) = 0. (3.20) 

We consider (|3.20|) as an operator equation on the Banach space X given by (|3.19|) . and we apply 
Theorem 13.51 with d 2 as the bifurcation parameter. Set v := u — u. Then (|3.20j) assumes the 
equivalent form 

Dn xx + An + G(v) = 0, v x (0)=v x (l) = 0. (3.21) 
where A is the Jacobian matrix of F evaluated at u and 



G(v) = F(v + u) - Av, G(0) = 0,G v (0) 



0. 



(3.22) 



Denote the left hand side of (|3.21[) by T(v, d 2 ), where T is a one-parameter family of operators acting 
on X and taking its elements into Y = C°([0, /];M 2 ). Clearly, T is a C 2 mapping. The spectrum of 
the linear operator Lio = T v (0, d 2c «t) = ^^(0, cfocrit) consists of the eigenvalues Hij of the matrices 
Bj given by (|3.9p with its corresponding eigenfunctions are ipj(x)yij, i = 1, 2, j = 0, 1, 2, • • • , where 
ifjj is given by (|3.7|) and y^- is the eigenvector of the matrix Bj corresponding to the eigenvalues /ijj 
(see (|3.8p ). Now, all matrices Bj = A — QD are to be taken at cfo = dicrit- As it can be seen from 



the proof of Theorem 13.51 and from (|3.10p for i = 1,2; for all nonnegative integer j except j = 1, 
all Hij have negative real parts. For j = 1 one eigenvalue pin is equal to and the other /i2i is 
negative. The eigenfunction corresponding to fin = is vi = yn cos(7ra/7). Thus, the null-space 
of the operator L±q = T v (0,d2 Cr it) is a one-dimensional linear space spanned by vi. Owing to the 
orthogonality and completeness of the eigenfunction system of the operator —-§^z, the range of this 
operator is given by 



R(L 10 ) = {w € C°([0,/];1R 2 ) : the eigenfunction expansion of 

TTX i r 7TX i 

w does not contain cos — j U span |y2i cos — | , 



so that the codimension of R(Liq) is one. 
Let L 12 = ^(0,d 2crii ). Then 



/ d 2 

L\2 = D — -k where D 
ox 



dD 
dd~2~ 




1 



Clearly, 



£l2Vl 



2 TTX i 

cos — D yn 



\ 2 TTX 





) cos T 




. V2 _ 



Under the assumption L12V1 jf( y 2 i cos ^p, we see that L12V1 does not belong to R(L\q), fulfilling 
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the condition (hi) of Theorem 13.61 
Letting 

Z = R(L\o), 

which is a closed subspace of Y, we verify that all the hypotheses of Theorem 13.61 hold; moreover, 
(0, d2crit) is a bifurcation point, and there exist a 5 > 0, a function d 2 : (—5,5) — ► M such that for 
s G (-5,5) 

v(x; s) = syn cos — + s0(x; s) 

is a solution of (|3.2ip with d% = ^(s), |s| < 5, 6^(0) = O,0(x;O) = 0, and d 2 G C 1 ,(^(x;-) G 
C 1 , 0(-; s) G Z. □ 

Remark 3.9. T/ie corresponding solution of 13.20\) . i.e. the non-constant stationary solution of 
the nonlinear parabolic system U.10\) is 

u(x; s) = u + syn cos ™ + 0(s 2 ), (3.23) 

(corresponding to the choice d 2 = d 2 (s), \ s\ < 5), i.e. 

N(x) = 7V + s??icos^ + 0(s 2 ), (3.24a) 
P(x) = P + sr] 2 cos ^ + 0(s 2 ). (3.24b) 

since s is considered to be small here, this solution is called as a small amplitude pattern. 

Remark 3.10. Because of Theorem \3.6\ $1. 10\) has no other stationary solution apart from (N,P) 
and (|3.23p in a neighborhood of (u, d 2cr it) Glxl. 

Remark 3.11. In the linear case (by Theorem \3. 5\) for the function d 2 holds: d 2 (s) = d 2cr u, an d 
a corresponding one parameter family of solutions is u + svi, s£R. 



4. Numerical approximation 

4-1. The numerical scheme 

The reaction-diffusion equations (jl.lOp are solved numerically using the forward Euler method 
in time, the centered difference method in space. This numerical scheme gives a stable solution 
under a certain that stasisfies the CFL (Courant-Friedrichs-Lewy) condition. The details are as 
follows. 

Consider the computational domain [0, 1] and the mesh size h and the time step size At, which 
will be determined later in (|4.10p . Set = p Denote by Nj and Pj the numerical approximation 
of N(jh, kAt), P(jh, kAt), respectively for j = 0, 1, ■ • ■ , and k = 1, 2, ■ • • . Then, given initial 
data N®, Pj,j = 0, 1, • • ■ , Nh, the numerical scheme is to solve 



n; 



k+l 



Nj + AtNj 



Pf + AtePf 



N k 

3 



aP k 

3 



pk 

3 



+ Nl 



+ Aid 



N}_ 



2N? 



+ 



N k , 



h 2 



7 + 5(3Pf 

TTppf 



+ 



N} 



pk + N k 



pk 

+ Atd 2 ^ 



2Pf 



+ ^1 



h 2 



(4.1a) 
(4.1b) 
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for j = 1, 2, ■ • ■ , Nh — 1, iteratively for fc = 1,2,--- . On the boundaries x = 0, x = 1 where Neumann 
condition holds, we used a three-point interpolation scheme to guarantee the second-order accuracy 
in space as follows: 

iVf _ 4JV* + 3j\r fc = 0; P 2 fc - 4Pf + 3P fc = 0; (4.2a) 

N k Nh - 2 ~ 4iVVi + = 0; P k Nh _ 2 - 4P^_! + 3P^ = 0. (4.2b) 

We will then establish the the positivity of the numerical solutions and boundedness for the numer- 
ical prey solution under certain conditions on At. Suppose that < Nj < 1 for j = 1, • • • , Nh — 1. 
Then, for j = 2, • • • ,N h -2, 



n; 



k+1 



Nf + AtNf 



< iVJ + AtiVJ 



1 - N. 



aP k 

3 



pk , N k 

3 3 



+ Atdi 



Nj_ x - 2iVj + N]f +1 



1-iV* 



+ 2^(1-^) 



Nf + At(l - iV*) 



AT* 2C? 1 

JV? H - 

< N k + At(l-N k )(l + ^) 



* , 2di 



/i 2 



< i-At(l + -^) + At(l + -^)<l 



(4.3) 



provided 1 — At(l + ^) > 0. For j = 1, owing to the boundary condition (|4.2|) . N k+l is given by 



fe+i 
l 



1-iVT 



Hence, the same analysis as above yields, instead of (|4.3|) . 



+ Atdi 



4(-iVf + iVf) 
3^2 ' 



(4.4) 



iV^ 1 < 



1 - At(1 + s?> 



JVf + At(l + < 1 



provided 1 — At(l + ^4-) > 0. Analgously, one gets 



iV^ii = AT^_i + AiiVVi 



1 AT fc aP ^-l 
1 - ^N h -1 ~ ~pk 



pk , pjk 



+ Atdi 



K-N k Nh ^+N k Nh _ 2 ) 
3/i 2 



and therefore 



1-A*(1 + ^)W_ 1 + A ( (1 + §)<1 



provided 1 - At(l + |^) > 0. On the other hand, suppose that < N k < 1 for j = !,-■■ ,N h -l. 
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Then, for j = 2, • • • , N h - 2, 



n; 



A'+l 



+ AtN* 



aP fc 

1 _ jyk 3 — 



H-Aidi^- 1 J J+ 



/; 2 



Aid 



> 2vJ + At ( 1 - ivj ) - Ai a JVf - 2 

> (1 + Ai - Aia)A^ - AiiVj 5 - 2^^A^ 



/i 2 J 



-, a a 2di 
1 - Aia - Ai— i 

n 2 



iV* > 0, 



(4.5) 



provided 1 — Aia — Ai 2 ^- > 0. Next for j = 1, by using (|4.4p . the procedure to get the estimate 
1.51) leads to 



1 - Aia - Ai 



Adi 
3^2 



iVf > 



(4.6) 



provided 1 — Aia — Ai|^ > 0. Similarly, under the same conditions, one obtains 

Mi' 



N k+1 > 



1 - Aia - Ai 



3/i 2 



N 



N h -1 



> 0. 



(4.7) 



Next, suppose that Pf > and < iVj < 1 for j = 1, ■ ■ ■ , N h - 1. Recalling (TOT) , one then 
obtains, for j = 2, • • • , iV/j — 2, 



P 



k+l 



P/ + AiePf 



6 — j 
S+- ^tt + 



i + ppf Pf + 



+ Aid 



> (1 - eoAi - ^?-)Pf > 



(4.8a) 



provided 1 — eoAi — 2 ^ d2 > 0. For j = 1 and j = iV/j — 1, taking into account of the boundary 
condition (|4.2p . one gets 



P*+i > (i _ e ^Ai - > forj = 1 and j = N h - 1 



(4.9a) 



provided 1 — eoAi — 4 ^f 2 > 0. Collecting all the above results, we are now in a position to state 
the following theorem: 



Theorem 4.1. Let < Nf < 1, < P? /or j = 0, • • • , iV^. Suppose that 



At < min 



/; 2 



/ ? 2 



a/i 2 + 2di ' /i 2 + 2di ' e<5/i 2 + 2d 2 



(4.10) 
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Then the numerical solutions Nj and Pj obtained iteratively by (|4.ip and (|4.2p satisfies that 

0<N}<1, 0<P/, /orj = l,--- ,JV h -l, /or A: = 0,1,2,.-. . (4.11) 

Numerically a steady state is declared to reach when either the L 2 or L max -noTva difference is 
less than a given tolerance value. The Li and L max -norm differences are defined as follows: 

\\u(-,kAt)\\l = / \{u steady (x,kAt) - u h (x,kAt)\ 2 dx, 
Jo 

||u(-, kAt)^ = max \u steady (x, kAt) - u h (x,kAt)\ , 

x6[0,i] 

where VL s teady are given by (|3.24|) with 0(s 2 ) terms neglected and Uh(x, kAt) is the piecewise linear 
interpolation of the numerical solution {Nj, Pj),j = 0, • • • , Nh- 

4-2. Numerical examples 

Set e = 1, a = 1.1,7 = 0.05,/? = 1,5 = 0.5. The unique positive equilibrium is (N,P) = 
(0.113585, 0.471397). If we fix I = 1 for the length of the habitat the interval (j3TTj) becomes 

1.488790091 x 10~ 3 < di < 5.95160365 x 10~ 3 . 

In the following Figure [H stability regions, the mean prey-predator diffusion coefficients, d\ and 
d2, are plotted. 

We tested our model in the cases of (di,d 2 ) = (0.005,0.2) and (di,d 2 ) = (0.005,0.32), which 
are in the stable and unstable regions with varying s = 0.05, 0.1, 0.2, 0.3, 0.4, respectively. In these 
cases, the critical value for Turing bifurcation d c is 0.271. Figure [2] shows the numerical prey and 
predator solutions, N and P, with respect to time at a specified fixed point x = 0.25. As shown in 
Figure[2j for (di,d 2 ) =(0.005,0.2), the equilibrium solution (N,P) is asymptotically stable and for 
(di, d 2 ) = (0.005, 0.32), the equilibrium solution (N, P) is unstable. For the simulation in the case of 
(dx, d 2 ) = (0.005, 0.2), we used the spatial mesh size h = 0.005, and the time step size At = 0.00006 
determined by the (|4.10p . The iteration was run until the time equals to 1000, with approximately 
1.6 • 10 iterations. In the case of (di,d 2 ) = (0.005,0.32), the mesh size /i=0.005 and the time 
step size At= 0.0000375 were used, which were alsothe (|4.10p . In this case also the simulation was 
done until the time equals to 1000, with approximately 2.6 - 10 7 iterations. In Figure [3l in case 
of (di,d 2 ) = (0.005,0.2), the prey and predator solutions are plotted with respect to number of 
iterations and space. We clearly see that as time goes to infinity, the solution converges to the 
equilibrium solution (N,P). In the lower figure in Figure El in case of (di,d 2 ) = (0.005,0.32), 
where d 2 is in unstable region, the prey and predator solutions are plotted with respect to number 
of iterations and space. We clearly see that as time goes to infinity, the solution shows the deviation 
from the equilibrium solution (iV, P). 

In Figured for the values near d c , (di,d 2 ) = (0.005,0.27) and (di,d 2 ) = (0.005,0.272) are 
considered. By varying s values from 0.05 to 0.4, the prey predator solution has a small amplitude 
pattern which we expected in the theory. In Figure [5] and Figure [61 we have plotted the prey and 
predator solutions and their small amplitude patterns with respect to number of iterations and 
space by changing s values. Near the d c , in case of (di,d 2 ) = (0.005,0.27), we use the mesh sizes 
h = 0.005, At = 0.0000444 and ran our simulation until the number of iteration is approximately 
10 7 . In case of (di,d 2 ) = (0.005,0.272), we have used with the mesh sizes h = 0.005 and At = 



14 



0.0000441176. Again our runs were continued until the number of iteration was approximately 10 . 
In Figure [5] and Figure [6l the axis scale in s = 0.1 has been used as that of the case of s = 0.4 
which has a bigger amplitude pattern. Comparing the solutions in Figure [5] and Figure [6] with the 
non-constant stationary solution (13. 24ft . we clearly observe that as time goes to infinity the prey 
and predator solutions converge to non-constant stationary solution (|3.24[) which confirms that 
(TV, P) undergoes a Turing bifurcation. 

Discussions 

System (jl.9p describes the dynamics of a ratio-dependent predator-prey interaction with diffu- 
sion. Prey quantity grows logistically in the absence of predation, predator mortality is neither a 
constant nor an unbounded function, but it is increasing with the predator abundance and both 
species are subject to Fickian diffusion in a one-dimensional spatial habitat from which and into 
which there is no migration. It is assumed that the system without diffusion has a positive equi- 
librium and under certain conditions it is asymptotically stable. We show that analytically at a 
certain critical value a diffusion driven (Turing type) instability occurs, i.e. the stationary solution 
stays stable with respect to the kinetic system (the system without diffusion). We also show that 
the stationary solution becomes unstable with respect to the system with diffusion and that Turing 
bifurcation takes place: a spatially non-homogenous (non-constant) solution (structure or pattern) 
arises. A first order approximation of this pattern (|3.23p is explicitly given. A numerical scheme 
that preserve the positivity of the numerical solutions and the boundedness of prey solution is 
introduced. Numerical examples are also included. 
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d 1 versus d 2 




0.002 0.003 0.004 0.005 



di 

Figure 1: d\ and di plot, from equation (|3.18p 
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numerical prey(N) solution 



numerical predator(P) solution 




Figure 2: Left: The prey solution at x—0.25 with respect to time, the constant line represents N e = N and the 
two solid lines represent two different d,2 values. Right: The predator solution at £=0.25 with respect to time, the 
constant line represents P e = P and the two solid lines represent two different d,2 values. 
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prey(N) 's time solution 



predator(P)'s time solution 



prey pattern N(x,t) predator pattern P(xJ) 




prey(N) 's time solution predator(P) 's time solution 

prey pattern N(x,t) predator pattern P(xJ) 




Figure 3: (di :0. 005,^2:0. 2, d c :0. 271) Upperleft: The prey solution N(x, t) with respect to time and space when d,2 < d c . 
Prey pattern shows the convergence to the equilibrium solution iV as time increases. Upperright: The predator 
solution P(x,t) with respect to space when d,2 < d c . Predator pattern shows the convergence to the equilibrium 
solution P as time increases. (di :0. 005, di :0. 32, d c :0. 271) LowerLeft: The prey solution N(x,t) with respect to time 
and space when d-2 > d c . Prey pattern shows the deviation from the equilibrium solution N as time increases. 
LowerRight: The predator solution P(x,t) with respect to space when di > d c . Predator pattern shows the deviation 
from the equilibrium solution P as time increases. 
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prey pattern(d=0.27) 



predator pattern(d=0.27) 
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prey pattern(d-0.272) 



predator pattern(d=0.272) 
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I 
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s=0.3 
s=0.4 



0.2 0.4 0.6 0.8 1 

space 



Figure 4: Upper: The prey/predator solution pattern N(x,t), P(x,t) when d,2 < d c with varing 
s. (di:0. 005,^2:0. 27, d c :0. 271) Lower: The prey/predator solution pattern N(x,t), P(x,t) when di > d c . 
(di:0.005,d 2 :0.272,d c :0.271) 
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predator(P)'s time solution 

predator pattern P(x,t) 









— ~tt2 
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predator(P)'s time solution 
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0.472 
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0.476 

0.475 

ii4'4 
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Figure 5: Upper Left: The predator solution pattern P(x,t) when s — 0.1, < d c . Upper Right: The predator 
solution pattern P(x,t) when s — 0.4, di < d c . (di:0. 005,^2:0. 27, d c :0. 271) Lower Left: The predator solution pat- 
tern P(x,t) when s = 0.1,^2 > d c . Upper Right: The predator solution pattern P(x,t) when s = 0.4, d^ > d c . 
(di:0.005,d 2 :0.272,d c :0.271) 
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prey(N)'s time solution 



prey(N)'s time solution 





Figure 6: Upper Left: The prey solution pattern N(x,t) when s — 0.1,^2 < d c . Upper Right: The prey solution 
pattern N(x, t) when s = 0.4, d 2 < d c . (di:0. 005,^2:0. 27, d c :0. 271) Lower Left: The prey solution pattern N(x, t) when 
s = 0.1, d 2 > d c . Upper Right: The prey solution pattern N(x,t) when s = 0.4, d 2 > d c . (di:0.005,d 2 :0.272,d c :0.271) 
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Figure 1 d\ and d 2 plot, from equation (|3.18|) 



Figure 2 Left: The prey solution at x=0.25 with respect to time, the constant line represents 
N e = N and the two solid lines represent two different d 2 values. Right: The predator solution 
at x=0.25 with respect to time, the constant line represents P e = P and the two solid lines 
represent two different d 2 values. 

Figure 3(e?i:0.005,<i2:0.2,<i c :0.271) Upperleft: The prey solution N(x,t) with respect to time 
and space when d 2 < d c . Prey pattern shows the convergence to the equilibrium solution 
N as time increases. Upperright: The predator solution P(x, t) with respect to space when 
d2 < d c . Predator pattern shows the convergence to the equilibrium solution P as time 
increases. {d\:Q. 005, d 2 :{). 32, d c :0. 271) LowerLeft: The prey solution N(x,t) with respect to time 
and space when d 2 > d c . Prey pattern shows the deviation from the equilibrium solution N 
as time increases. LowerRight: The predator solution P(x,t) with respect to space when 
c?2 > d c . Predator pattern shows the deviation from the equilibrium solution P as time 
increases. 

Figure 4 Left: The prey/predator solution pattern N(x,t), P(x,t) when d 2 < d c with varing 
s. (c?i:0.005,(i2:0.27,d c :0.271) Right: The prey/predator solution pattern N(x, t), P(x, t) when 
d 2 > dc. (di:0.005,d 2 :0.272,d c :0.271) 

Figure 5 Upper Left: The predator solution pattern P(x,t) when s = 0.1, d 2 < d c . Upper 
Right: The predator solution pattern P(x, t) when s = 0.4, di < d c . (di:0.005,d2-'0-27,d c :0.271) 
Lower Left: The predator solution pattern P(x,t) when s = 0.1, d 2 > d c . Upper Right: The 
predator solution pattern P(x,t) when s = 0.4,^2 > d c . (di:0.005,<i2:0.272,(i c :0.271) 

Figure 6 Upper Left: The prey solution pattern N(x,t) when s = 0.1,^2 < d c . Upper Right: 
The prey solution pattern N(x,t) when s = 0.4, d 2 < d c . (di:0.005,<i2:0.27,(f c :0.271) Lower 
Left: The prey solution pattern N(x,t) when s = 0.1,^2 > d c . Upper Right: The prey 
solution pattern N(x,t) when s = 0A,d 2 > d c . (di:0.005,d 2 :0.272,d c :0.271) 
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